Method and System for Adaptive Maximum Intensity Projection Ray Casting

ABSTRACT

The adaptive MIP ray casting system first fragments a 3-D dataset into multiple sub-volumes and constructs an octree data structure with each sub-volume being associated with one node of the octree data structure. The system then establishes a 2-D image plane and selectively launches a plurality of rays towards the 3-D dataset, each ray adaptively interacting with a subset of the sub-volumes and identifies the maximum data value along the ray path. The maximum data value is then converted into a pixel value on the 2-D image plane. Finally, the system interpolates pixel values at those locations where no pixel value is generated by ray casting and thereby generates a 2-D image of the 3-D dataset.

CROSS-REFERENCE RELATED APPLICATIONS

The present application is a continuation of application serial no.10/569,325, filed Sep. 15, 2006, which application is a U.S. NationalPhase Application of International Application PCT/US04/26819, filedAug. 18, 2004, which application is a continuation of application Ser.No. 10/642,797, filed Aug. 18, 2003, all of which applications arehereby incorporated herein by reference in their entirety.

FIELD OF THE INVENTION

The present invention relates generally to the field of volume datarendering, and more particularly, to a method and system for rendering avolume dataset using adaptive data representation, ray casting, and datainterpolation techniques.

BACKGROUND OF THE INVENTION

Visualization of the internal structure of 3-D objects on a 2-D image isan important topic within the field of computer graphics and has beenapplied to many industries, including medicine, geoscience,manufacturing, and drug discovery.

For example, a CT scanner can produce hundreds or even thousands ofparallel 2-D image slices of a patient's body including differentorgans, e.g., a heart, each slice including a 2-D array of data valuesand each data value representing a scalar attribute of the body at aparticular location, e.g., density. All the slices are stacked togetherto form an image volume or a volumetric dataset of the patient's bodywith the heart embedded therein. A 2-D image showing the 3-D structuralcharacteristics of the heart is an important aid in the diagnosis ofcardiovascular disease.

As another example, the oil industry uses seismic imaging techniques togenerate a 3-D image volume of a 3-D region in the earth. Some importantgeological structures, such as faults or salt domes, may be embeddedwithin the region and not necessarily on the surface of the region.Similarly, a 2-D image that fully reveals the 3-D characteristics ofthese structures is critical in increasing oil production.

Maximum intensity projection (MIP) ray casting is a technique developedfor visualizing the interior of a solid region represented by a 3-Dimage volume on a 2-D image plane, e.g., a computer monitor. Typically,a plurality of rays are cast from a 2-D radiation plane into the imagevolume, each ray casting responsible for identifying the maximum datavalue (e.g., intensity) at a voxel within the image volume along itsrespective ray path and transferring it into an image value at a pixelon a 2-D image plane through a predefined screen transfer function. Theimage value is indicative of the 3-D characteristics of the objectsembedded within the image volume encountered by the ray path, e.g.,their shapes and orientations. The image values associated with thepixels on the 2-D image plane form an image that can be rendered on acomputer monitor.

Going back to the CT example discussed above, even though a doctor canarbitrarily generate 2-D image slices of the heart by intercepting theimage volume in any direction, no single image slice is able tovisualize the whole surface of the heart. In contrast, a 2-D imagegenerated through MIP ray casting of the CT image volume can easilyreveal the 3-D characteristics of the heart, which is very important inmany types of cardiovascular disease diagnosis. Similarly in oilexploration, MIP ray casting of 3-D seismic data can help petroleumengineers to determine more accurately the 3-D characteristics ofunderground geological structures of a region that are potential oilreservoirs and to increase oil production significantly.

Even though MIP ray casting plays a key role in many important fields,there are several technical challenges that need to be overcome toassure wide deployment of the MIP ray casting technique. First, MIP raycasting is a computationally expensive process. In order to produce ahigh quality 2-D image that can capture the 3-D characteristics of a 3-Dtarget, MIP ray casting needs to process a large 3-D dataset, whichusually means a large number of calculations. For example, it requiresat least 140 million calculations to generate a 2-D image of 512² pixelsfor a typical 3-D dataset of 512³ voxels using conventional MIP raycasting algorithms.

Moreover, many applications require that MIP ray casting of a 3-Ddataset operate in real-time so that a user is able to view successive2-D images of the 3-D dataset, each 2-D image having different viewingangles or visualization parameters, without a significant delay. Inmedical imaging, it is generally accepted that a sequential 2-D imagerendering of at least six frames per second meets the need for real-timeinteractive feedback. This is equivalent to nearly 1 billioncalculations per second.

Given the limited computational capacity of modern computers, moreefficient algorithms have been developed to reduce the computationalcost of MIP ray casting. However, many of these algorithms achieve theirefficiency by sacrificing the quality of the generated 2-D images. Forexample, a common problem with discrete representation of a continuousobject is the jitter effect, which is most obvious when a user zooms into view more details of the continuous object. If the jitter effect isnot carefully controlled, it may significantly corrupt the quality of animage generated by a MIP ray casting algorithm.

Therefore, it would be desirable to develop a new MIP ray casting methodand system that increase the rendering efficiency while having less orpreferably imperceptible impact on the image quality.

SUMMARY OF THE INVENTION

A preferred embodiment of the present invention is an adaptive MIP raycasting method and system that generate high-quality 2-D images of 3-Dobjects embedded in a 3-D dataset using adaptive data representation,MIP ray casting, and data interpolation techniques.

The method and system first fragment the 3-D dataset into multiplesub-volumes of different sizes and construct a data structure, e.g., anoctree, associating nodes of the octree with different sub-volumes ofthe dataset, each sub-volume also having a set of data value parameterscharacterizing the data value distribution within the sub-volume. Inparticular, the whole dataset is associated with the root node of theoctree.

The method and system then cast a plurality of rays from a 2-D radiationplane towards the 3-D dataset. In one embodiment, commonly referred toas parallel projection, the plurality of rays are parallel to one otherand each ray has a unique ray origin on the 2-D radiation plane. Inanother embodiment, commonly referred to as perspective projection, theplurality of rays are all launched from the same ray origin and each rayhas a unique launching angle with respect to the 3-D dataset.

Each ray, having a certain cross-section or ray thickness and apredefined current data value record, interacts with a series ofsub-volumes of the 3-D dataset along its ray path and updates thecurrent data value record whenever it identifies a higher data valuealong the ray path. Specifically, during each interaction between theray and a sub-volume, one or more comparisons are performed between thecurrent data value record and a set of data value parameters associatedwith the sub-volume. If the current data value record is smaller thanthe data value at a voxel of the sub-volume on the ray path, e.g., themiddle point of the ray path within the sub-volume, the data value atthe voxel will replace the old record and become the current data valuerecord. This process continues until the ray exits the 3-D dataset orcertain predefined ray casting termination conditions have beensatisfied, identifying the maximum data value record of the 3-D dataseton the ray path. The maximum data value record is translated into apixel value on a 2-D image plane in accordance with a predefined screentransfer function. The pixel values on the 2-D image plane associatedwith the plurality of rays constitute a 2-D image illustrating the 3-Dinternal structures embedded in the 3-D dataset.

There are two competing requirements for the choice of an appropriatesub-volume that interacts with a ray. On the one hand, a largesub-volume can significantly reduce the computational cost of the raycasting, rendering the method and system more efficient; on the otherhand, a small sub-volume can characterize the 3-D internal structuresmore accurately, producing an image of higher quality. The method andsystem has to strike a balance between these two competing goals byselecting a sub-volume that meet certain predefined conditions, e.g.,that the difference between the maximum and minimum data values of asub-volume is smaller than a predefined threshold. In one embodiment,the whole dataset associated with the root node of the octree is firstselected to interact with the ray. However, since the set of data valueparameters associated with the root node rarely captures the 3-Dinternal structures embedded in the 3-D dataset sufficiently, thedataset often has to be further fragmented into smaller sub-volumes, andone of the smaller sub-volumes that is closest to the ray origin alongthe ray path and associated with a child node of the root node is thenexamined to determine if it meets the predefined conditions. Thisfragmentation and examination process repeats recursively until asub-volume of desired size is identified.

Once a sub-volume is finally chosen for interacting with a ray, themethod and system estimate the data value at a particular location inthe sub-volume that is on the ray path through tri-linear interpolationof known data values surrounding that particular location. If theestimated data value is greater than the ray's current data valuerecord, the estimated data value is assigned to the ray as its newcurrent data value record. Otherwise, the ray's current data valuerecord remains unchanged. Following that, if the ray is still inside the3-D dataset, the method and system identify another sub-volume ofappropriate size next to the chosen one and farther along the ray path;a new estimated maximum data value is determined and the estimationcontinues.

The computational cost and the image resolution of a MIP ray castingalgorithm are determined primarily by the number and density of rayscast from a 2-D radiation plane. Given a limited number of rays that canbe handled by modern computer hardware within a reasonable time period,the method and system optimize the locations of the ray origins on the2-D radiation plane such that the pixel values associated with thoserays on the 2-D image plane can effectively capture the characteristicsof the 3-D internal structures embedded within the 3-D dataset.

In one embodiment, four rays are first launched from four corners of the2-D radiation plane. After getting the four pixel values at the fourcorners of the 2-D image plane, the method and system check thevariation of the values against a predefined imaging error threshold. Ifthe variation is above the threshold, one more ray is launched from thecenter of the radiation plane, thereby dividing the 2-D image plane intofour sub-planes. Variations of the pixel values within each sub-planeare further checked against the predefined imaging error threshold. Anysub-plane whose variation is above the threshold is further dividedrecursively until the variation of any sub-plane is below the threshold.

Finally, the method and system estimate a value at each pixel on the 2-Dimage plane that is not associated with a ray based on the surroundingpixel values on the 2-D image plane. In one embodiment, those pixelvalues are results of bi-linear interpolation of pixel values at thefour corners of a smallest sub-plane that have met the predefinedimaging error threshold.

BRIEF DESCRIPTION OF THE DRAWINGS

The aforementioned features and advantages of the invention as well asadditional features and advantages thereof will be more clearlyunderstood hereinafter as a result of a detailed description ofpreferred embodiments of the invention when taken in conjunction withthe drawings.

FIGS. 1A and 1B are two schematic illustrations of two embodiments ofthe present invention, parallel projection and perspective projection,respectively.

FIG. 1C depicts three 2-D images of a 3-D object at different imageresolutions.

FIG. 1D depicts a 3-D object sampled by a 3-D dataset.

FIG. 2 is a flowchart depicting a MIP ray casting algorithm according toone embodiment of the present invention.

FIG. 3A depicts one embodiment of a 3-D dataset fragmentation.

FIG. 3B depicts a sub-volume comprising eight cells or voxels and scalardata value distribution within a voxel or a cell.

FIG. 4 is a flowchart illustrating a computer implementation of the 3-Ddataset fragmentation.

FIG. 5 schematically illustrates a series of interactions between a rayand sub-volumes of the 3-D dataset.

FIG. 6 is a flowchart depicting a computer program for identifyingmaximum data value along a ray path according to an embodiment of thepresent invention.

FIG. 7 presents an 1-D example illustrating the adaptive MIP ray castingaccording to one embodiment of the present invention.

FIG. 8 is a flowchart illustrating how a computer program adaptivelygenerates a 2-D image according to some embodiments of the presentinvention.

FIG. 9 presents an example of bi-linear interpolation of pixel values.

FIGS. 10A and 10B illustrate two methods of adjusting a predefinedimaging error threshold.

FIG. 11 is a block diagram of an adaptive MIP ray casting systemaccording to one embodiment of the present invention.

FIG. 12 is a block diagram of the adaptive MIP ray casting renderingsystem implemented a computer cluster according to one embodiment of thepresent invention.

FIG. 13 shows two images of a human body through different versions ofadaptive MIP ray casting algorithm according to one embodiment of thepresent invention.

DESCRIPTION OF EMBODIMENTS

FIGS. 1A and 1B schematically illustrate two embodiments of the presentinvention, parallel projection and perspective projection, respectively.In a 3-D domain represented by the Cartesian coordinates (x, y, z),there is an image volume 106 containing one or more 3-D objects. MIP raycasting generates a 2-D image by casting a plurality of rays 104 intothe image volume 106 and identifying the maximum data value along eachray path so as to visualize the 3-D characteristics of the 3-D objectson the 2-D image. Note that the shape of a ray in the context of thepresent invention is not a 1-D line, but a 3-D tube or cone depending onspecific ray configurations discussed below.

In the embodiment shown in FIG. 1A, the plurality of rays 104 arelaunched from different locations on a radiation plane 102, each raytraveling in parallel with other rays 104 towards the image volume 106.The shape of a ray is a tube 103 that has a cross-section of constantsize, also referred to as ray thickness. Different rays correspond todifferent pixels on a 2-D image plane 101, and each pixel value isdetermined by the maximum data value of the image volume identified on acorresponding ray path. Such MIP ray casting configuration is called“parallel projection”. In the case of parallel projection, the imageplane 101 can be located anywhere along the ray path. Illustratively,the image plane 101 can be co-incident with the radiation plane 102 asshown in FIG. 1A.

In the embodiment shown in FIG. 1B, the plurality of rays 104 arelaunched from an identical ray origin 107 on the radiation plane 102,each ray having a unique perspective angle with respect to the imagevolume 106. As a result, the shape of each individual ray is a cone 105that has a cross-section of variable size that is a function of both aviewing angle 108 and a distance 109 between the ray origin 107 and theimage volume 106. Similarly, different rays correspond to differentpixels on a 2-D image plane 101, and each pixel value is determined bythe maximum data value of the image volume identified on a correspondingray path. Such MIP ray casting configuration is called “perspectiveprojection”. In the case of perspective projection, the image plane 101can be located anywhere along the ray path between the ray origin 107and the image volume 106. Illustratively, the image plane 101 can beco-incident with the radiation plane 102 as shown in FIG. 1B.

Generally speaking, the image resolution of a 2-D image generated by oneof the two MIP ray casting configurations discussed above depends on thenumber of pixels per unit area. In practice, the total number of pixelsper image is often a constant, e.g., 512², which is dependent upon thenumber of pixels on a computer monitor. Therefore, higher imageresolution in parallel projection is achieved by reducing thecross-section of a ray tube 103 and thereby increasing the number ofpixels per unit area, and moving the radiation plane 102 along the raypaths away from or towards the image volume 106 does not affect theimage resolution. In contrast, higher resolution in perspectiveprojection is achieved by either reducing the viewing angle 108 or thedistance 109 or both.

FIG. 1(C) illustrates how the image resolution changes as a function ofray thickness. For simplicity, each image has 10²=100 pixels, and theray thicknesses of grids 110-1, 112-1 and grid 114-1 are S, 4S and 16S,respectively. As a result, the image generated on grid 114-2 has thelowest resolution, while the image generated on grid 110-2 has thehighest resolution; and the image generated on grid 112-2 hasintermediate resolution. Parallel projection achieves such magnificationeffect by reducing the size of the radiation plane and therebyincreasing the density of ray origins on the radiation plane; andperspective projection achieves the same effect by moving the cameracloser to the 3-D objects in the image volume or changing the camera'sperspective angle.

The image volume 106 is usually represented by a 3-D dataset 116 asshown in FIG. 1D. Each voxel of the 3-D dataset has a scalar data valueindicative of a physical attribute of the 3-D objects at a particularlocation. As shown in FIG. 1D, the 3-D dataset 116 has three orthogonalcoordinates x, y and z. Along each coordinate, the dataset comprises astack of 2-D data slices perpendicular to that coordinate. Each dataslice includes a 2-D array of data values on a regular grid along theother two coordinates. For example, the three orthogonal data slices116-1, 116-2 and 116-3 provide different (but limited) perspectives ofimage volume 106 in the directions z, x, and y, respectively.

Algorithm Overview

MIP ray casting according to one embodiment of the present inventioncomprises three major steps as shown in FIG. 2. At step 201, the methodfragments an image volume into multiple sub-volumes of different sizes,each sub-volume being associated with a set of data parameterscharacterizing the data value distribution of the sub-volume. Detailsabout this method are described below in connection with FIGS. 3A-3B andFIG. 4.

At step 203, the method performs a plurality of ray castings against theimage volume from a radiation plane that is located a certain distanceaway from the center of the image volume and oriented in a certaindirection with respect to the image volume. A ray cast from theradiation plane selectively interacts with a series of sub-volumeslocated along its ray path. The maximum data value on the ray path isconverted into a pixel value in the image plane. More details about thisprocess are provided below in connection with FIGS. 5-7.

At step 205, the method estimates pixel values at any other locations inthe image plane using the pixel values generated at step 203 andpresents a 2-D view of the 3-D objects embedded in the image volume.Steps 203 and 205 may be repeated numerous times to estimate more pixelvalues. In addition, steps 203 and 205 may be repeated to generate aseries of new images when the radiation plane's location or directionchanges or different visualization parameters are applied. Furtherdetails about this process are provided below in conjunction with FIGS.8-10.

Image Volume Fragmentation

According to one embodiment of the present invention, an image volume isfirst fragmented into a set of sub-volumes and each sub-volume isfurther fragmented into smaller sub-volumes. Such recursivefragmentation continues until the smallest sub-volumes reach predefinedsize limits. All the sub-volumes, including the whole image volume, areassociated with a hierarchical data structure, which provides a newrepresentation of the original image volume.

FIG. 3A depicts a 3-D dataset 302 after being fragmented into aplurality of sub-volumes (e.g., 585 sub-volumes). These sub-volumesbelong to four fragmentation levels I, II, III, and IV. For illustrativepurposes, sub-volumes of different sizes, e.g., II-2, III-2 and IV-2,have been removed from FIG. 3A to expose other sub-volumes of smallersizes. Each of the sub-volumes is associated with one node of an octreedata structure 304 based on its size. Note that only a small subset ofthe sub-volumes are shown in the octree 304 to simplify theillustration.

Starting from fragmentation level I, the whole dataset 302 is treated asa single sub-volume I-1 and associated with the root node of octree 304.At fragmentation level II, dataset 302 is partitioned into eightsub-volumes of the same size; and each sub-volume, e.g., II-6, isassociated with one intermediate node at level II of octree 304. Atfragmentation level III, each sub-volume such as sub-volume II-6 isfurther divided into eight smaller sub-volumes III-1, . . ., and III-8;and similarly, each sub-volume, e.g., III-6, is associated with anintermediate node at level III. At fragmentation level IV, eachsub-volume, such as sub-volume III-6, is partitioned into eightsub-volumes (including sub-volume IV-6), and each sub-volume at level IVis associated with a leaf node of octree 304. Since there are 8⁰sub-volumes at level I, 8¹ sub-volumes at level II, 8² sub-volumes atlevel III, and 8³ sub-volumes at level IV, the dataset 302 isfragmented, into a total of

8⁰+8¹+8²+8³=1+8+64+512=585

sub-volumes at four different fragmentation levels.

A sub-volume has an associated set of data parameters characterizing thedata value distribution within the sub-volume. In one embodiment, theset of parameters includes at least three elements:

-   -   V_(min) representing the minimum of the data value within the        sub-volume;    -   V_(avg) representing the average of the data value within the        sub-volume; and    -   V_(max) representing the maximum of the data value within the        sub-volume.

As illustrated in FIG. 3A, the image volume fragmentation is a recursiveprocedure. This procedure does not terminate until the smallestsub-volumes reach a predefined size limit. In one embodiment, thispredefined size limit or the smallest sub-volume associated with anoctree leaf node is a sub-volume having 2×2×2 cells or voxels of theimage volume. A cell or voxel is the smallest unit of the image volumethat has an associated data value. For example, if the image volume has512³ cells which is typical in MIP ray casting and the smallestsub-volume has 2×2×2 cells, the image volume fragmentation processcontinues to fragmentation level IX and the number of smallestsub-volumes is 256³.

FIG. 3B shows a sub-volume 306 having eight cells or voxels, two cellsor voxels next to each other in each of the three orthogonal directions.A voxel 308 is a cube that is centered at one sampled data value, e.g.,V₀. The data value within a voxel is assumed to be constant and theremay be a discontinuity of data value across the voxel's boundary. Incontrast, a cell 310 is a cube of the same size that has eight sampleddata values V₁-V₈, one at each corner of the cube. According to oneembodiment, the data value at any point of the cell is approximated bytri-linear interpolation and there is no data value discontinuity acrossthe cell's boundary. Note that the present invention applies equally toimage volumes represented in either format, cell or voxel. Forsimplicity, the following discussion focuses on the cell representationas an example.

Referring again to FIG. 3B, there is a 3-D Cartesian coordinate systemwith its origin at the center of cell 310 and the lengths of the cell'sthree orthogonal edges are D_(x), D_(y), and D_(z), respectively. Thedata value at any point of the cell 310, V(x,y,z), can be expressed bytri-linear interpolation as a function of the data values at the eightcorners of the cell 310, V₁-V₈, as follows:

${V\left( {x,y,z} \right)} = {{\frac{V_{1}}{8}\left( {1 + \frac{2\; x}{D_{x}}} \right)\left( {1 + \frac{2\; y}{D_{y}}} \right)\left( {1 + \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{2}}{8}\left( {1 + \frac{2\; x}{D_{x}}} \right)\left( {1 - \frac{2\; y}{D_{y}}} \right)\left( {1 + \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{3}}{8}\left( {1 - \frac{2\; x}{D_{x}}} \right)\left( {1 - \frac{2\; y}{D_{y}}} \right)\left( {1 + \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{4}}{8}\left( {1 - \frac{2\; x}{D_{x}}} \right)\left( {1 + \frac{2\; y}{D_{y}}} \right)\left( {1 + \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{5}}{8}\left( {1 + \frac{2\; x}{D_{x}}} \right)\left( {1 + \frac{2\; y}{D_{y}}} \right)\left( {1 - \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{6}}{8}\left( {1 + \frac{2\; x}{D_{x}}} \right)\left( {1 - \frac{2\; y}{D_{y}}} \right)\left( {1 - \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{7}}{8}\left( {1 - \frac{2\; x}{D_{x}}} \right)\left( {1 - \frac{2\; y}{D_{y}}} \right)\left( {1 - \frac{2\; z}{D_{z}}} \right)} + {\frac{V_{8}}{8}\left( {1 - \frac{2\; x}{D_{x}}} \right)\left( {1 + \frac{2\; y}{D_{y}}} \right){\left( {1 - \frac{2\; z}{D_{z}}} \right).}}}$

FIG. 4 is a flowchart summarizing a computer program for image volumefragmentation as well as octree construction discussed above. At step401, the computer program receives a 3-D dataset and initializes anoctree data structure by associating the 3-D dataset with the root nodeof the octree. At step 403, the dataset is fragmented into eightsub-volumes and eight new child nodes are generated at the root nodesuch that each sub-volume is associated with one of the eight childnodes.

Starting from step 405, each sub-volume is recursively fragmented intosmaller sub-volumes until the smallest sub-volumes reach a predefinedsize limit. At step 405, the computer program first checks if the sizeof a sub-volume reaches the predefined size limit, e.g., a sub-volume of2×2×2 cells. If not, at step 407, the sub-volume is fragmented intoeight smaller sub-volumes. One of the eight smaller sub-volumes ispicked and the recursive fragmentation starts at step 405 again.

If the size of a sub-volume has reached the predefined size limit, e.g.,2×2×2 cells, the computer program stops fragmenting this sub-volume atstep 409 and the node of the octree associated with the sub-volume isone of the leaf nodes of the octree. At step 411, the computer programchecks if there are other sub-volumes that have not been fragmented to2×2×2 cells. If true, one of those sub-volumes is selected and therecursive fragmentation continues at step 405. Otherwise, the imagevolume fragmentation is terminated.

As discussed above, each sub-volume is associated with a set of dataparameters such as (V_(min), V_(avg), V_(max)) for estimating a pixelvalue on the 2-D radiation plane. In one embodiment, after the 3-Ddataset has been fragmented and the octree has been constructed, thecomputer starts from the leaf node level of the octree, calculating thecorresponding (V_(min), V_(avg), V_(max)) of a sub-volume associatedwith a leaf node, such as sub-volume IV-6 in FIG. 3A. Since a sub-volumeassociated with a leaf node comprises only eight cells and each cell haseight data values, most of which are shared by two or more cells, thecomputer program only needs to retrieve 27 data values from a storagedevice, e.g., a hard disk, to calculate (V_(min), V_(avg), V_(max)) ofthe sub-volume.

After processing all the sub-volumes at the leaf node level, thecomputer program moves one level up on octree 304 to calculate thecorresponding parameters for a sub-volume like III-6 that is associatedwith an intermediate node. This process continues until it reaches theroot node. Since each intermediate node has eight child nodes and eachchild node is associated with a sub-volume whose associated dataparameters have been determined previously, each data value parameter ofthis intermediate node can be expressed as a function of the eightparameters associated with the eight child nodes:

-   -   V_(min)=Min(V_(min) _(—) ₁, V_(min) _(—) ₂, . . . , V_(min) _(—)        ₈);    -   V_(avg)=(V_(avg) _(—) ₁+V_(avg) _(—) ₂+ . . . +V_(avg) _(—)        ₈)/8; and    -   V_(max)=Max(V_(max) _(—) ₁, V_(max) _(—) ₂, . . . , V_(max) _(—)        ₈).

As will be apparent, this generation of data value parameters proceedsin a direction opposite to the octree construction which starts at theroot node level in a top-down fashion as shown in FIG. 3A. Thisbottom-up approach of data parameters generation is most efficient sinceit maximally reuses the calculation results at a lower level of theoctree.

After step 201, a new representation of the original image volume isavailable for step 203 of the MIP ray casting. This representationincludes an octree and a plurality of sub-volumes of different sizes,each sub-volume associated with one node of the octree and having a setof data parameters characterizing the data value distribution within thesub-volume.

3-D Adaptive MIP Ray Casting

MIP ray casting is an algorithm that generates a 2-D image revealinginternal 3-D structures embedded in an image volume, and each pixelvalue of the 2-D image is estimated by casting a ray from the pixel intothe image volume, identifying the maximum data value along the ray pathand converting the maximum data value into a pixel value in accordancewith a predefined screen transfer function. According to one embodimentof the present invention, an efficient 3-D adaptive MIP ray castingalgorithm includes three steps:

-   -   Step one: the algorithm constructs a mathematical model for a        ray having a predefined current data value record and a certain        cross-section, and represents a 3-D object as a discrete 3-D        dataset, each data value of the dataset representing a physical        property of interest at a particular location of the object.    -   Step two: the algorithm compares a series of data values along        the ray path with the ray's current data value record, one at a        time, and replaces the current data value record if a higher        data value is identified.    -   Step three: the algorithm converts the ray's final data value        record, which is also the highest data value along the ray path,        into a pixel value on a 2-D image in accordance with a        predefined screen transfer function.

The screen transfer function is typically a monotonic function that ischosen arbitrarily to illuminate different aspects of the internalstructure embedded in a 3-D dataset. FIG. 5 is a screen transferfunction that may be used in MIP ray casting. The vertical axisrepresents the magnitude of pixel values on a 2-D image and thehorizontal axis represents the magnitude of data values of a 3-Ddataset. The screen transfer function can be divided into threesections: any data values smaller than V_(low) correspond to a constantpixel value P_(low), any data value greater than V_(high) correspond toa constant pixel value P_(high), and any data values between V_(low) andV_(high), e.g., V_(i), correspond to a pixel value P_(i) that is betweenP_(low) and P_(high).

FIG. 6 is a flowchart illustrating how a computer program identifies thehighest data value of an image volume along a particular ray pathaccording to one embodiment of the present invention. Since the imagevolume has been fragmented into a plurality of sub-volumes that areassociated with an octree as shown in FIG. 3A, each sub-volume having aset of data parameters (V_(min), V_(avg), V_(max)) characterizing thedata value distribution therein, the computer program begins theidentification process with the whole image volume associated with theroot node of the octree for efficiency purposes.

At step 601, the computer program selects a ray from a plurality of rayslaunched from a 2-D radiation plane, the ray having an origin on the 2-Dradiation plane, a specific orientation with respect to the imagevolume, and an initial current data value record V_(rec). In oneembodiment, the initial current data value record is set to be V_(low)of the screen transfer function as shown in FIG. 5, since any data valuelower than V_(low) corresponds to the same pixel value P_(low) on a 2-Dimage. As discussed below in connection with FIG. 7, this approach, whenapplied to the octree data structure, can effectively eliminate asignificant amount of the image volume from the identification process,making the identification process more efficient.

At step 603, the computer program selects from the octree data structurethe largest sub-volume that encounters the ray for the first time tocheck if the ray path within this sub-volume includes a point having adata value higher than the ray's current data value record V_(rec). Bydefault, the first sub-volume selected at step 603 is always the wholeimage volume, which has an associated set of data parameters (V_(min),V_(avg), V_(max)). However, as shown in FIG. 6, the computer program mayre-visit step 603 after completing steps 608 or 619 if some conditionshave been met. In this case, the computer program should have examinedat least one sub-volume associated with the octree data structure, andthe largest sub-volume selected at step 603 will not be the whole imagevolume but the sub-volume that is spatially next to a sub-volume fromwhich the ray has just exited.

As mentioned above, the goal of the computer program is to identify thehighest data value on the ray path within the image volume. However, toreduce imaging error as well as to expedite the identification process,the ray's current data value record is not replaced by a higher valuejust because the sub-volume selected at step 603 has at least one datavalue higher than V_(rec), e.g., V_(rec)<V_(max). Rather, a sub-volumeneeds to meet three criteria to replace the ray's current data valuerecord:

-   -   the maximum data value V_(max) of the sub-volume must be higher        than the ray's current data value record V_(rec) (step 605);    -   the maximum screen value differential (MSVD) within the        sub-volume must be below a predefined threshold (step 609); and    -   the data value V_(est) at a predefined location on the ray path        within the sub-volume, e.g., the middle point of the ray path        within the sub-volume, must be higher than the ray's current        data value record V_(rec) (step 615).

At step 605, the computer program compares the ray's current data valuerecord V with the maximum data value V of the sub-volume selected atstep 603. If V_(rec)>V_(max), there is no chance for the ray's currentdata value record to be updated. The computer program accordinglyeliminates this sub-volume and the ray keeps its current data valuerecord when exiting the sub-volume. Note that the ray may have traveledthrough the whole image volume just after exiting the sub-volume.Therefore, the computer program checks if this is the last sub-volume onthe ray path at step 608. If so, the computer program stops the raycasting and the ray's current data value record V_(rec) will be treatedas the highest data value along the ray path, which is subsequently usedfor estimating a corresponding pixel value on the 2-D radiation plane.Otherwise, if V_(rec)<V_(max), there is at least one location in thesub-volume whose value is higher than the ray's current data valuerecord. But it requires some further steps taken by the computer programto determine if this voxel can be used to update V.

At step 607, the computer program estimates the maximum screen valuedifferential (MSVD) of the sub-volume using the predefined screentransfer function (STF). The MSVD of a sub-volume is indicative of theimpact of the data value variation within the sub-volume on the pixelvalue variation on the 2-D radiation plane. In one embodiment, dependingupon the magnitudes of the sub-volume's associated data parametersV_(min), V_(avg), V_(max) and the current data value record V_(rec), theMSVD of the sub-volume may have slightly different definitions:

-   -   MSVD=MAX(STF(V_(max))−STF(V_(avg)), STF(V_(avg))−STF(V_(min)))        if V_(rec)<V_(min);    -   MSVD=MAX(STF(V_(max))−STF(V_(avg)), STF(V_(avg))−STF(V_(rec)))        if V_(min)<V_(rec)<V_(avg); or    -   MSVD=STF(V_(max))−STF(V_(rec)) if V_(avg)<V_(rec)<V_(max).

Since the purpose of the identification process is to identify thehighest data value on the ray path within an image volume, any voxel ofa sub-volume whose value is below the current data value record V_(rec)should have no bearing on the estimation of the MSVD of the sub-volume.Meanwhile, V_(avg) is more accurate than V_(max) and V_(min) in terms ofcharacterizing the overall data value distribution within thesub-volume, and the MSVD definition that involves V_(avg) as well asV_(max) and V_(min) can produce more accurate images than the MSVDdefinition that only involves V_(max) and V_(min). In some alternativeembodiments, the average data value V_(avg) may be replaced by a mediandata value V_(med). The median data value is the one that is higher thanone half of the sub-volume's data values and lower than another half ofthe sub-volume's data values.

At step 609, the computer program compares the estimated MSVD with thepredefined threshold. If the estimated MSVD is above the threshold, thecomputer program steps down along the octree data structure from thenode associated with the current sub-volume and selects a smallersub-volume that is associated with one of the child nodes (step 611).Among the multiple smaller sub-volumes associated with the child nodes,the computer program selects the one that encounters the ray ahead ofall the other sub-volumes and then returns to step 605 to determine ifthe newly selected sub-volume is small enough to pass the test at step609. Note that the newly selected sub-volume has its own set of dataparameters (V′_(min), V′_(avg), V′_(max)) and, in particular, itsmaximum data value V′_(max) may be smaller than the maximum data valueV_(max) of the parent node. If V′_(max)<V′_(rec), the computer programwill skip this smaller sub-volume and move to step 608 as discussedabove.

Following the loop including steps 605-611, the computer program mayreach a leaf node of the octree data structure that is associated with asmallest sub-volume having, e.g., 2×2×2 cells. If the smallestsub-volume still fails to pass the test at step 609, the computerprogram may select one of the eight cells or even fragment the cell intoeight sub-cells and then examine one of the sub-cells in order to passthe test at step 609. As shown in

FIG. 3B, any data value at any point within a cell or sub-cell can beestimated using the eight data values at the eight corners through,e.g., tri-linear interpolation. On the other hand, a ray, according tosome embodiments, is deemed to be a tube or a cone having a certaincross-section. Therefore, the fragmentation within the cell has to stopat a certain stage, e.g., if the size of a sub-cell is smaller than thecross-section of the ray within the sub-cell. If this happens, thecomputer program may assume the test at step 609 has been met and moveforward to step 613 automatically.

Assuming that the test at step 609 has been met or the MSVD of asub-volume (or a cell or sub-cell) is below the predefined threshold,the computer program then compares the current data value record V_(rec)with an estimated data value V_(est) at one point on the ray path withinthe sub-volume to determine whether or not V_(rec) should be updated byV_(est). In some embodiments, the center point of the ray path withinthe sub-volume is chosen for this comparison; and in some otherembodiments, the computer may choose positions such as the locationwhere the ray enters/exits the sub-volume. Since the chosen position isoften not exactly one corner of a cell, the computer program often hasto interpolate the data value V_(est) at a specified position. To do so,the computer program at step 613 first determines which cell containsthe chosen position and then retrieves the data values at the eightcorners of the cell from the computer memory or hard disk to estimatethe data value V_(est) at that position through, e.g., tri-linearinterpolation.

If V_(rec)>V_(est) (“No” at step 615), the ray keeps its current datavalue record while exiting the current sub-volume. Note that wheneverthe ray exits a sub-volume, the computer program needs to check at step608 if the ray has also exited the whole image volume. IfV_(rec)<V_(est) (“Yes” at step 615), the computer program replaces thecurrent data value record V_(rec) with V_(est) at step 617 and V_(est)becomes the ray's current data value record. Next, the computer programat step 619 compares the ray's current data value record with V_(high)of the screen transfer function, because, as shown in FIG. 5, a datavalue record above V_(high) does not have any impact on thecorresponding pixel value P_(high). If V_(rec)>V_(high), there is noneed for continuing the ray casting and the computer program terminatesthe ray casting immediately for this particular ray. Otherwise, the raycasting will continue. However, since the ray has just exited asub-volume, the computer program needs to determine if the ray isoutside the whole image volume first before resuming ray casting at step603.

The adaptive MIP ray casting algorithm shown in FIG. 6 is essentially aprocedure of selecting a sequence of sub-volumes along a ray path andidentifying one position in one of the sub-volumes that has the highestdata value along the ray path. Since different sub-volumes areassociated with different nodes of an octree data structure andcharacterized by different sets of data parameters (V_(min), V_(avg),V_(max)), the algorithm can adaptively eliminate those sub-volumes thatare irrelevant to the goal of the MIP ray casting (e.g., at step 605)and significantly reduce the computational cost.

FIG. 7 presents an 1-D example illustrating how the algorithm works inidentifying a voxel having the highest value among a group ofarbitrarily distributed data values according to the embodiment shown inFIG. 6. Note that the goal in 1-D is easier to accomplish than that in2-D or 3-D, because the 2-D or 3-D version of the algorithm needs toestimate a sequence of data values on different sections of a ray pathwithin different sub-volumes and identify the highest value among them.In contrast, the highest data value is a known to the 1-D version of thealgorithm and what is unknown is the location of this highest value,while the location of the highest value on the ray path is not a concernof the 2-D or 3-D version of the algorithm (even though it is readilyavailable).

In this example, there are sixteen data samples along the horizontalaxis, A through P. For example, sample A is 100, sample J is 221 andsample N is 85, etc. An intuitive approach is to first assign thecurrent data value record an arbitrary value (e.g., 180) and thencompare it with all the sixteen data samples one by one. Whenever asample has a higher value, the current data value record will bereplaced by the higher value. This process continues until the lastsample P is examined and the sample that has the highest value isidentified to be sample I (227).

Apparently, the computational cost of this intuitive approach is O(N), Nbeing the number of data samples. If there are M² groups of data samples(analogous to M² ray castings required for generating a 2-D image), thetotal computational cost will be O(NM²). If N is a very small number,this approach might work. However, in a real-life application, e.g., CTscan image volume, N is often a large number (e.g., a few hundred oreven more than a thousand), rendering O(NM²) too high to be handled bymost computer systems.

Instead, a 5-level binary tree 700 is constructed such that each of thesixteen data samples is associated with one of the sixteen leaf nodes ofthe binary tree 700. Every non-leaf node of the binary tree 700,including the root node 740, is associated with a set of data parameters(V_(min), V_(avg), V_(max)) indicative of the minimum, average andmaximum data values at the leaf nodes that are descendants of thisnon-leaf node. For example, node 720 has a set of data parameters (100,112, 125), since it has only two leaf nodes whose associated values are100 and 125, respectively, and the root node 740 has a set of dataparameters (58, 138, 227), because all the sixteen leaf nodes are itsdescendants and among the sixteen leaf nodes, leaf node I has thehighest value 227 and leaf nodes L and M have the lowest value 58.

Assuming that the initial data value record V_(rec)=180, V_(low)=180,V_(high)=220, respectively, the predefined threshold T_(MSVD)=10 and thescreen transfer function STF(V)=180 if V<V_(low) or V ifV_(low)<V<V_(high) or 220 if V_(high)<V. The root node is first selectedto compare with V_(rec) according to step 603. Since V_(rec) (180) issmaller than V_(max) (227) at the root node, there is at least one ofthe sixteen samples having an associated value higher than 180.Accordingly, the MSVD at the root node is estimated according to step607 as (STF(V_(max))−STF(V_(rec)))=227−180=47, because V_(rec) (180) islarger than V_(avg) (138) at the root node. Since MSVD (47) is greaterthan T_(MSVD) (10), the computer program has to investigate a sub-groupof leaf nodes by stepping one level down to one of the child nodesaccording to step 611. Since V_(max) of the child node 760 is only 171which is smaller than the current data value record V_(rec) (180), thereis no need to consider any of the eight leaf nodes that are descendantsof the child node 760 according to step 605. Instead, the child node 780is chosen since it is next to the child node 760 and its V_(max) (227)is greater than V_(rec) (180). This process continues until the node 784is identified which satisfies all the three criteria discussed above andthe value (227) at leaf node I is chosen to replace the current datavalue record (180). Note that even though there are still six leaf nodes(K to P) that are not examined, there is no need for any furtherinvestigation because the current data value record V_(rec) (227) ishigher than V_(high) (220) and any higher value will have no impact onthe resulting pixel value. Therefore, the computer program mayappropriately terminate this process according to step 619.

Mathematically, the computational cost of the binary tree approachdiscussed above is only O(log₂N) if the data samples are randomlydistributed, which is substantially lower than the computer cost of theintuitive approach O(N). When the algorithm is extended to 3-D byfragmenting a 3-D dataset into multiple sub-volumes and associating eachsub-volume with a node of an octree data structure as discussed above,the computational cost can drop substantially.

The predefined threshold used by the adaptive MIP ray casting at step609 is a parameter affecting the efficiency and accuracy of thealgorithm. Given a predefined screen transfer function, the thresholddictates how many sub-volumes the computer program needs to examine atstep 609 until it finds the sub-volume that can pass the test. In otherwords, the larger the threshold, the more efficient the algorithm.Meanwhile, the algorithm at step 613 estimates the data value at anarbitrarily selected position, e.g., the center point of the ray pathwithin a sub-volume, treats this estimated value as the maximum datavalue on the ray path within the sub-volume and compares it with thecurrent data value record at step 615. This arbitrarily selectedposition, however, may not be the exact position at which the data valuereaches its peak on the ray path within the sub-volume. A smallerthreshold translates into a smaller sub-volume that can pass the test atstep 609 and therefore a smaller error between the arbitrarily selectedposition and the exact position. In practice, the computer program hasto strike a balance between these two competing goals, efficiency versusaccuracy, by selecting an appropriate threshold.

According to one embodiment, the threshold used at step 609 is not astatic value, but a dynamic one that varies from one image to anotherimage, to generate an optimal result at a reasonable cost. Assuming thatE₀ is an initial threshold defined by a user, a customized threshold isdefined as:

-   -   E(E₀, P_(zoom))=E₀/Sqrt(P_(zoom)) in the case of a parallel        projection, where P_(zoom) is the zoom factor indicating the        physical size of an image; and    -   E(E₀, P_(distance))=E₀*Log₂(P_(distance)) in the case of a        perspective projection, where P_(distance) is the distance        between the ray origin and the center of the image volume.

Increasing the zoom factor P_(zoom) in the case of parallel projectiondecreases the threshold, and the computer program needs to test moresub-volumes before identifying the maximum data value along a ray path,generating images with a higher resolution and less jitter effect. Incontrast, increasing the distance P_(distance) between the ray originand the image volume in the case of perspective projection has theopposite impact on the image quality.

Besides image quality, user-friendliness is another factor whenevaluating an adaptive MIP ray casting algorithm. For example, any timeinterval between the time the radiation plane is positioned with respectto the image volume and the image is rendered on the computer monitor isless desirable in dynamic imaging applications like cardiovasculardisease diagnosis. To reduce such time delay and offer the user morecontrol during the diagnosis process, in one embodiment, the predefinedthreshold can be further adjusted by the ratio between a desired frameper second (DFPS) and an actual frame per second (AFPS):

-   -   E(E₀, P_(zoom), DFPS, AFPS)=E(E₀, P_(zoom))*DFPS/AFPS in the        case of a parallel projection; and    -   E(E₀, P_(distance), DFPS, AFPS)=E(E₀, P_(distance))*DFPS/AFPS in        the case of a perspective projection.

As indicated by their names, DFPS means an image rendering speedpreferred by a user while AFPS is the one that is limited by the imagingsystem's settings such as the computer hardware, the number of pixels onthe radiation plane and the image volume being processed. If thepreferred image rendering speed is higher than the one offered by thesystem, the threshold should increase accordingly. As a result, anincreased threshold achieves the user-preferred image rendering speed bycompromising the image quality. On the other hand, if the preferredimage rendering speed is lower than the one that the system can provide,the system may have additional capacity for lowering the threshold andgenerating higher-quality images at a user-acceptable speed.

2-D Image Estimation

A 3-D adaptive MIP ray casting identifies the maximum data value along aray path in an image volume and then converts it into a pixel value at aparticular location on a 2-D image plane according to a predefinedscreen transfer function. Since the 3-D MIP adaptive ray casting is acomputationally expensive operation, the efficiency of the algorithm issubstantially dependent on the number of pixel values on the image planethat have to be generated through ray casting. FIG. 8 is a flowchartillustrating how a computer program adaptively determines which pixelvalues on the image plane need to be generated through the adaptive MIPray casting process and which pixel values can be estimated using thosegenerated pixel values with negligible, if not imperceptible, effects onthe image quality.

At step 801, the computer program defines a 2-D image plane that is acertain distance away from and oriented in a certain direction withrespect to an image volume as shown in FIGS. 1A and 1B. The image planeis usually represented by a 2-D array in the computer memory includingM×M (e.g., 512×512) elements, each element of the 2-D array storing onepixel value at a corresponding location on the image plane.

At step 803, the computer program subdivides the image plane into aplurality of mini-planes of the same size, each mini-plane correspondingto a portion of the 2-D array and including P×P (e.g., 16×16) pixels. Inone embodiment, a mini-plane is associated with the root node of aquadtree data structure.

At step 805, for each mini-plane created at step 803, the computerprogram conducts four separate adaptive MIP ray castings and estimatesfour pixel values at the four corners of the mini-plane.

Starting from step 807, the computer program recursively fills everyelement of the 2-D array with a pixel value such that, at the end ofthis recursive process, a 2-D image of the image volume is formed in the2-D array. In one embodiment, the computer program adaptively subdividesa mini-plane into multiple smaller sub-planes, each sub-plane associatedwith one child node of the quadtree data structure and assigns eachpixel of the sub-plane a pixel value through adaptive MIP ray casting oran interpolation method. More specifically, the computer program firstchecks if there is any mini-plane that has not been processed at step807, i.e., if there is any mini-plane comprising at least one pixelwithout a pixel value. If not, every pixel of the image plane hasalready been assigned a pixel value and the computer program terminates.Otherwise, the computer program examines a quadtree data structureassociated with a mini-plane identified at step 807 for any sub-planethat has not been processed at step 809.

If all the sub-planes associated with the quadtree data structure havebeen processed, the computer program returns to step 807 to process thenext mini-plane. If not, the computer program at step 811 identifies asub-plane and checks if this sub-plane comprises only 2×2 pixels. Asub-plane having only 2×2 pixels is the smallest sub-plane that amini-plane can be subdivided into. In one embodiment, the computerprogram estimates a pixel value for any pixel on the sub-planeidentified at step 811 through adaptive MIP ray casting at step 813. Ifthe sub-plane identified at step 811 is larger than a sub-plane of 2×2pixels, the computer program needs to conduct further experiments todecide how to fill the sub-plane with pixel values.

At step 815, the computer program first estimates the maximum pixelvalue differential (MPVD) of this sub-plane. As indicated by its name,the MPVD measures the pixel value distribution within a sub-plane. Inone embodiment, the MPVD is defined as

MPVD=MAX(|P _(avg) −P ₁ |,|P _(avg) −P ₂ |,|P _(avg) −P ₃ |,|P _(avg) −P₄|)/S,

where S is the total number of pixels within a sub-plane, P₁-P₄ arepixel values at the corners of the sub-plane and P_(avg) is the averageof P₁-P₄

$P_{avg} = {\frac{1}{4}{\left( {P_{1} + P_{2} + P_{3} + P_{4}} \right).}}$

In some other embodiments, the average pixel value P_(avg) may bereplaced by the median pixel value P_(med) of the sub-plane for thepurpose of estimating the MPVD. The MPVD is compared with a predefinedimaging error threshold for determining how to fill the remaining pixelsin the sub-plane that have no associated pixel values.

Generally speaking, interpolation is a more efficient alternative thanadaptive MIP ray casting and therefore should be employed wheneverpossible. On the other hand, a pixel value generated through adaptiveMIP ray casting is more accurate than that generated throughinterpolation. Since different portions of the 2-D image capturedifferent 3-D structures in the image volume and may require differentlevels of imaging resolution, a universal imaging error threshold acrossthe 2-D image is not the most favorable option.

Instead, at step 817, the computer program selectively adjusts thethreshold when dealing with different portions of the image volume. Forexample, the computer program should choose a smaller predefined imagingerror threshold when rendering those portions of the 2-D image that maycorrespond to the edge of certain 3-D internal structure, e.g., a bonein medical imaging, which subsequently forces a sub-plane to be furtherfragmented into smaller sub-planes and, as a result, more pixel valuesin the sub-plane have to be generated through adaptive MIP ray castinginstead of interpolation. More details about the customization of thepredefined imaging error threshold modification are provided below inconjunction with FIG. 10.

At step 819, the computer program checks if the MPVD of the sub-plane isabove the customized imaging error threshold. If not, the computerprogram first estimates a pixel value at the center of the sub-planethrough adaptive MIP ray casting and then subdivides the sub-plane intomultiple smaller sub-planes at step 823, each smaller sub-planeassociated with one child node of the quadtree data structure. If theMPVD is higher than the customized imaging error threshold, the computerprogram assigns to any unfilled pixel a value through interpolation ofexisting pixel values at step 821. In one embodiment, the pixel valuedistribution within a sub-plane is approximated by a bi-linearinterpolation of the four pixel values at the corners of the sub-plane.As shown in FIG. 9, the pixel value P(x,y) at any location (x,y) can bebi-linearly interpolated as

${P\left( {x,y} \right)} = {{\frac{P_{1}}{4}\left( {1 + \frac{2\; x}{d_{x}}} \right)\left( {1 + \frac{2\; y}{d_{y}}} \right)} + {\frac{P_{2}}{4}\left( {1 - \frac{2\; x}{d_{x}}} \right)\left( {1 + \frac{2\; y}{d_{y}}} \right)} + {\frac{P_{3}}{4}\left( {1 - \frac{2\; x}{d_{x}}} \right)\left( {1 - \frac{2\; y}{d_{y}}} \right)} + {\frac{P_{4}}{4}\left( {1 + \frac{2\; x}{d_{x}}} \right){\left( {1 - \frac{2\; y}{d_{y}}} \right).}}}$

After step 821 or 823, the computer program returns to step 809 toprocess the next sub-plane, if any. This process continues until everyelement of the 2-D array is allocated a pixel value corresponding to onelocation on the image plane that is generated through either 3-Dadaptive ray casting or bi-linear interpolation. As a result, a 2-Dimage of the image volume is formed that can be rendered on a graphicaldisplay device such as a computer monitor.

In general, the human eye is more sensitive to the high-frequencycomponents of an object, e.g., a 2-D image. However, bi-linearinterpolation (FIG. 9) tends to smear out the high frequency componentsof an image because it is essentially a low-pass filter. Therefore, tocapture the high-frequency components and control the totalcomputational cost, the computer program may generate more pixel valuesthrough the more expensive adaptive MIP ray casting at certain portionswhile estimating most of the pixel values through bi-linearinterpolation at other portions of the image plane.

As discussed above, the predefined imaging error threshold can beadjusted to control the image resolution and image rendering speed. Forexample, if a sub-plane does not correspond to complex internalstructures or high-frequency components such as the edge of an object,the predefined imaging error threshold is raised to a high value andvice versa. However, to achieve this purpose, the computer program needsto predict whether a sub-plane being processed includes high-frequencycomponents or not.

According to some embodiments, two methods can be used to adjust thepredefined imaging error threshold to strike a balance between imagequality and image rendering speed. As shown at step 817 of FIG. 8, theyare:

-   -   bi-linear interpolation error analysis; and    -   octree traversal edge detection.

As shown in FIG. 10A, a sub-plane 1000 has four known pixel values P₁-P₄at its four corners that are generated through adaptive MIP ray casting.Therefore, the bi-linearly interpolated pixel value P_(5′) at the centerof the sub-plane 1000 is

$P_{5^{\prime}} = {\frac{1}{4}{\left( {P_{1} + P_{2} + P_{3} + P_{4}} \right).}}$

Meanwhile, the computer program also generates a pixel value P₅ throughadaptive MIP ray casting at the center of the sub-plane 1000 when it issubdivided into four smaller sub-planes. The absolute difference betweenthe two pixel values |P₅−P₅| shows how accurate the result of bi-linearinterpolation is within this sub-plane 1000. For example, a higherabsolute difference not only means the sub-plane needs to be furthersubdivided, but also suggests that there may be high-frequencycomponents within the sub-plane. The computer program, accordingly, maysignificantly lower the imaging error threshold to fully capture the 3-Dcharacteristics of the internal structures within the sub-plane. On theother hand, a smaller absolute difference suggests that the image withinsub-plane is smooth and the predefined imaging error threshold may beincreased accordingly to reduce the computational cost withoutsubstantially degrading the image resolution.

Octree traversal edge detection is a method of adjusting the imagingerror threshold by examining the size difference between a sub-volume ora cell or a sub-cell that has been identified as including the maximumdata value and an image plane to catch a smaller 3-D object or the sharpedge of a 3-D object embedded within an image volume.

For simplicity, FIG. 10B depicts a 2-D dataset 1010 that comprises atleast two objects 1040 and 1060 and the 2-D dataset includes N×N pixels.Inside the dataset is an image line 1020 that comprises L pixels. Alongthe ray path 1030 are a series of sub-planes of different sizes thathave met at least one of the criteria as shown in FIG. 6 (steps 605, 609and 615, respectively). A zoom-in view 1080 depicts details of the 2-Ddataset surrounding a sub-plane in which the computer program identifiesthe maximum data value 1090 along the ray path 1030. Assuming that thefragmentation level of the sub-plane containing the maximum data valueis Q, a factor Z indicating the existence of small objects or sharpedges of an object near the ray path is defined as:

Z=L−N/2^(Q).

A higher Q (i.e., a smaller sub-plane) results in a higher Z factor,suggesting that there may be some small objects along the ray path 1030,and the predefined imaging error threshold should be lowered accordinglyso as not to miss those small objects or sharp edges along the ray path.On the other hand, a lower Q (i.e., a larger sub-plane) results in alower Z factor, indicating that it is less possible to identify smallobjects along the ray path, and the predefined imaging error thresholdcan be increased without losing significant details on the image line1020. Note that the Z factor is useful in adjusting the predefinedimaging error threshold when the ray exits the dataset without beingearly terminated inside the dataset.

Finally, the predefined imaging error threshold can also be adjustedaccording to the ratio of desirable frame per second and actual frameper second. If the former is higher than the latter, the predefinedimaging error threshold may increase so that more image frames can berendered than the system's default capability with a slight loss ofimage resolution; otherwise, the predefined imaging error threshold maydecrease if the former is lower than the latter.

Computer System

The adaptive MIP ray casting method discussed above can be implementedin a regular computer system and does not need any special hardwaresupport. FIG. 11 illustrates such a system 1100 in accordance with someembodiments of the present invention. An adaptive MIP ray casting system1100 typically comprises one or more central processing units (CPU's)1102, memory 1114, and one or more communication buses 1112 forinterconnecting the various components of system 1100. The adaptive MIPray casting system 1100 also includes a user interface 1104, including,for example, a display device 1106 for displaying 2-D images of a 3-Ddataset and a keyboard 1108 for receiving a user's image renderingrequirements. The system 1100 may optionally include a network or othercommunications interface 1110 for retrieving 3-D datasets from remotestorage devices or transmitting rendering results to remote clients.

Memory 1114 includes high-speed random access memory and may alsoinclude non-volatile memory, such as one or more magnetic disk storagedevices (not shown). The memory 1114 may also include mass storage thatis remotely located from the central processing unit(s) 1102. The memory1114 preferably stores:

-   -   an operating system 1116 that includes procedures for handling        various basic system services and for performing hardware        dependent tasks;    -   a network communication module 1118 that is used for connecting        system 1100 to various security devices or client computers (not        shown) and possibly to other servers or computers via one or        more communication networks (wired or wireless), such as the        Internet, other wide area networks, local area networks,        metropolitan area networks, and so on;    -   a system initialization module 1120 that initializes other        modules and data structures stored in memory 1114 required for        the appropriate operation of system 1100;    -   an adaptive volume rendering engine module 1122 that generates        2-D images of a 3-D dataset and renders it on display device        1106;    -   a 3-D dataset 1140 and a corresponding Octree structure 1142        representing the fragmented 3-D dataset;    -   a 2-D image dataset 1146 and a corresponding Quadtree structure        representing the fragmented 2-D image; and    -   a screen transfer function table 1144 that converts a data value        of the 3-D dataset into a pixel value of the 2-D image dataset.

The MIP ray casting engine 1122 includes executable procedures,sub-procedures, and other data structures supporting the imagegenerating process, such as:

-   -   a 3-D dataset fragmentation module 1124 that fragments a 3-D        dataset into multiple sub-volumes, estimates a set of data value        parameters for each sub-volume, constructs an octree data        structure and associates the sub-volumes with various nodes on        the octree data structure;    -   a 3-D adaptive MIP ray casting module 1126 that selects a series        of sub-volumes that interact with a ray based on a set of        predefined criteria and then identifies the maximum data value        along the ray path; and    -   a 2-D image rendering module 1128 that subdivides a 2-D image        plane into multiple mini-planes, each mini-plane being        associated with a quadtree data structure, and invokes the 3-D        adaptive ray casting module 1126 selectively to constructs a 2-D        image of the 3-D dataset.

As shown in FIG. 8, during the 2-D image construction, the system 1100first subdivides an image plane into a plurality of mini-planes and thenprocesses each mini-plane recursively, associating every pixel in themini-plane with a pixel value through bi-linear interpolation oradaptive MIP ray casting, if necessary. At the end, a single 2-D imageis formed in the memory 1114. Since the imaging result of one mini-planeis independent from that of another mini-plane, the adaptive MIP raycasting algorithm is actually a parallel algorithm that can beimplemented on a parallel computer or a cluster of computers.

FIG. 12 illustrates one embodiment of the present invention that employsa computer cluster 1200. The cluster 1200 comprises:

-   -   a data storage device 1210, e.g., a hard disk, that stores one        or more image volumes;    -   a plurality of computer servers 1230-1, 1230-2, 1230-3, 1nd        1230-4, each server preferably having one or more CPUs and its        own memory storing the direct volume rendering software;    -   a computer monitor 1250 for displaying the images generated by        the cluster;    -   a keyboard 1270 and a mouse 1290 for receiving the command and        rendering parameters from a user; and    -   a communication bus 1220 connecting the various components.

Illustratively, there are four servers 1230-1, 1230-, 1230-3, and 1230-4in the cluster 1200, and each server is responsible for generating aquarter of the final image, as indicated by the partial images 1240-1,1240-2, 1240-3, and 1240-4. Within each server, the job may be furtherpartitioned between different CPUs, if there is more than one CPU perserver. After all the partial images 1240-1 through -4 have beengenerated, one of the servers assembles the partial images into acomplete image and renders it on the monitor 1250.

EXAMPLES

FIGS. 13A and 13B provide an example of two images generated from a 3-DCT scan dataset using adaptive MIP ray casting according to oneembodiment the present invention. The two images are substantiallysimilar except that a pixel value in FIG. 13A is only dependent upon itsrespective maximum data value identified along the ray path according toa predefined screen transfer function, while a pixel value in FIG. 13Bis further weighted by the distance between the image plane and thelocation of the maximum data value in the image volume. As a result, thesame maximum data values identified at different parts of the imagevolume may cause different pixel values on the image of FIG. 13B,offering certain 3-D perspectives to the viewer. For example, the twoimage sections 1320 and 1340 of the femur bone in FIG. 13A havevirtually the same brightness, while the same two sections 1360 and 1380in FIG. 13B have different brightness. In particular, since section 1360is farther away from the viewer than section 1380, it looks slightlydimmer when compared with section 1380.

The foregoing description, for purpose of explanation, has been setforth with reference to specific embodiments. However, the illustrativediscussions above are not intended to be exhaustive or to limit theinvention to the precise forms disclosed. Many modifications andvariations are possible in view of the above teachings. The embodimentswere chosen and described in order to best explain the principles of theinvention and its practical applications, to thereby enable othersskilled in the art to best utilize the invention and various embodimentswith various modifications as are suited to the particular usecontemplated.

1. A method of adaptive maximum intensity projection ray casting,comprising: fragmenting a sampled 3-D dataset of a scalar field into aplurality of sub-volumes of different sizes, each sub-volume associatedwith a set of data value parameters characterizing the data valuedistribution of the scalar field within the sub-volume; defining ascreen transfer function that is dependent upon data values of thescalar field; selectively casting a plurality of rays towards thesampled dataset, each ray having an initial data value record whereinthe step of selectively casting is performed on a computer processor andcomprises, for each ray, adaptively selecting a subset of the pluralityof sub-volumes of different sizes that are located along a path of theray based on the data value distribution of the scalar field within thesub-volumes that are located along the path of the ray; identifying amaximum data value on the ray path that is within the selected subset;and converting the maximum data value into a pixel value in a 2-D imageplane according to the screen transfer function; and wherein the step ofselectively casting a plurality of rays towards the sampled datasetfurther comprises: subdividing the 2-D image plane into a plurality ofmini-planes; and for each of the plurality of mini-planes, estimatingfour pixel values at four corners of the mini-plane; and recursivelysubdividing the mini-plane into multiple sub-planes, and estimatingpixel values at other locations in the sub-planes by casting raystowards the sampled dataset until a maximum pixel value differential ofeach sub-plane is below a predefined imaging error threshold; and usingthe pixel values determined for the cast rays to estimate other pixelvalues at other locations in the 2-D image plane.
 2. The method of claim1, wherein the step of fragmenting a sampled 3-D dataset includes:fragmenting the sampled dataset into eight sub-volumes; and for eachsub-volume, recursively fragmenting it into eight smaller sub-volumesuntil the size of the smallest sub-volumes reaches a predefined sizelimit.
 3. The method of claim 2, wherein the predefined size limit is asub-volume comprising 2×2×2 3-D cells, a data value of the scalar fieldbeing associated with each corner of each cell.
 4. The method of claim3, wherein each cell has eight corners and the data value at a locationwithin the cell is tri-linearly interpolated using the data values atthe eight corners of the cell.
 5. The method of claim 1 wherein the setof data value parameters associated with a sub-volume include a maximum,an average, and a minimum data value of the scalar field within thesub-volume.
 6. The method of claim 1, further comprising: constructingan octree data structure having a root node, a plurality of intermediatenodes and a plurality of leaf nodes; associating the root node with thesampled 3-D dataset; associating each of the plurality of leaf nodeswith one of the smallest sub-volumes; and associating each of theplurality of intermediate nodes with one of the sub-volumes that islarger than the smallest sub-volume.
 7. (canceled)
 8. The method ofclaim 1, wherein the maximum pixel value differential of a sub-plane isdefined as the maximum deviation of pixel values at the four corners ofthe sub-plane from the average pixel value of the sub-plane.
 9. Themethod of claim 1, wherein the predefined imaging error threshold isadjusted by an image rendering speed provided by a user, a distance toan edge of an object embedded in the 3-D dataset, and a differencebetween a pixel value estimated from an adaptive MIP ray casting and apixel value estimated from a bi-linear interpolation.
 10. The method ofclaim 1 wherein the step of selecting a subset of the plurality ofsub-volumes that are located along the ray path includes: identifying alargest sub-volume that encounters the ray along the ray path and itsassociated set of data value parameters; estimating a maximum screenvalue differential of the sub-volume using the ray's current data valuerecord and the sub-volume's associated set of data value parameters; andrecursively selecting a smaller sub-volume contained in the sub-volumeuntil the estimated maximum screen value differential of the smallersub-volume is below a predefined threshold.
 11. (canceled)
 12. Themethod of claim 1, wherein the step of identifying the maximum datavalue on the ray path that is within the selected subset of theplurality of sub-volumes includes estimating a data value at apredefined position of a sub-volume; updating the ray's current datavalue record if the estimated data value is higher than the current datavalue record; and repeating said estimating and updating steps againstevery remaining sub-volume of the subset until the ray exits the sampleddataset.
 13. The method of claim 12, wherein the sub-volume is asmallest sub-volume comprising 2×2×2 cells, or a cell or a sub-cellwhose dimension is at least above the crosssection of the ray.
 14. Themethod of claim 1, wherein the step of using the pixel values determinedfor the cast rays to estimate other pixel values at other locationsincludes for each of the other locations, selecting four pixel valuessurrounding the location and associated with four cast rays; andbi-linearly interpolating a pixel value at the location using the fourpixel values.
 15. An adaptive maximum intensity projection (MIP) raycasting system, comprising: one or more central processing units forexecuting programs; a user interface for receiving a plurality of MIPray casting parameters; and an adaptive MIP ray casting engine moduleexecutable by the one or more central processing units, the modulecomprising: instructions for fragmenting a sampled 3-D dataset of ascalar field into a plurality of sub-volumes of different sizes, eachsub-volume associated with a set of data value parameters characterizingthe data value distribution of the scalar field within the subvolume;instructions for defining a screen transfer function that is dependentupon data values of the scalar field; instructions for selectivelycasting a plurality of rays towards the sampled dataset, each ray havingan initial current data value record; instructions for adaptivelyselecting a subset of the plurality of sub-volumes of different sizesthat are located on a ray path based on the data value distribution ofthe scalar field within the sub-volumes that are located along the pathof the ray; instructions for identifying a maximum data value on the raypath that is within the selected subset; instructions for converting themaximum data value into a pixel value in a 2-D image plane according tothe screen transfer function; and wherein the instructions forselectively casting a plurality of rays towards the sampled datasetfurther comprise instructions for: subdividing the 2-D image plane intoa plurality of mini-planes; and for each of the plurality ofmini-planes, estimating four pixel values at four corners of themini-planes; and recursively subdividing the mini-plane into multiplesub-planes, and estimating pixel values at other locations in thesub-planes by casting rays towards the sampled dataset until a maximumpixel value differential of each sub-plane is below a predefined imagingthreshold; and instructions for using the pixel values determined forthe cast rays to estimate other pixel values at other locations in the2-D image plane.
 16. The system of claim 15, wherein the instructionsfor fragmenting a sampled 3-D dataset include: fragmenting the 3-Ddataset into eight sub-volumes; and for each sub-volume, recursivelyfragmenting it into eight smaller sub- volumes until the size of thesmallest sub-volumes reaches a predefined size limit.
 17. The system ofclaim 16, wherein the predefined size limit is a sub-volume comprising2×2×2 3-D cells, a data value of the scalar field being associated witheach corner of each cell.
 18. The system of claim 17, wherein each cellhas eight corners and the data value at any location within the cell istri-linearly interpolated using the data values at the eight corners ofthe cell.
 19. The system of claim 15, wherein the set of data valueparameters include a maximum, an average, and a minimum data value ofthe scalar field within the sub-volume.
 20. The system of claim 15,further comprises: instructions for constructing an octree datastructure comprising a root node, a plurality of intermediate nodes anda plurality of leaf nodes; instructions for associating the root nodewith the 3-D dataset; instructions for associating each of the pluralityof leaf nodes with a smallest sub-volume from the plurality ofsub-volumes; and instructions for associating each of the plurality ofintermediate nodes with a sub-volume from the plurality of sub-volumesthat is larger than the smallest sub-volume. 21-53. (canceled)